CD33 CAR NK: D1

suppressPackageStartupMessages({
library(tidyverse)
library(scran)
library(scater)
library(patchwork)
library(scDblFinder)
library(cowplot)
library(pheatmap)
library(edgeR)
library(ggrepel)
})

Loading Data

data.dir <- "Files_already_generated_with_Sevenbridges/D1/"

data1 <- read.csv(paste(data.dir, "_1_Combined_Rhapsody19_Target1_RSEC_MolsPerCell.csv", sep = ""), skip = 7)
data1$Cell_Index <- paste(data1$Cell_Index, "_S1", sep = "")

data2 <- read.csv(paste(data.dir, "Combined_Rhapsody19_Target2_RSEC_MolsPerCell.csv", sep = ""), skip = 7)
data2$Cell_Index <- paste(data2$Cell_Index, "_S2", sep = "")

data <- rbind(data1, data2) 

The metadata contains Sample_Tag (nucleotide labeled antibody used to distinguish between samples), Sample_Name (patient id), multiplet status assigned by the SevenBridges pipeline

# batch 1
## table with cell information (sample barcode, donor, ...)
coldata1 <- read.csv(paste(data.dir, "_1_Rhapsody19_Target1_Sample_Tag_Calls.csv", sep = ""), skip = 7) %>%
  dplyr::mutate(invalid = (Sample_Name == "Multiplet" | Sample_Name == "Undetermined")) %>%
  dplyr::mutate(batch = "S1")
coldata1$Cell_Index <- paste(coldata1$Cell_Index, "_S1", sep = "")

# batch 2
## table with cell information (sample barcode, donor, ...)
coldata2 <- read.csv(paste(data.dir, "Rhapsody19_Target2_Sample_Tag_Calls.csv", sep = ""), skip = 7) %>%
  dplyr::mutate(invalid = (Sample_Name == "Multiplet" | Sample_Name == "Undetermined")) %>%
  dplyr::mutate(batch = "S2")
coldata2$Cell_Index <- paste(coldata2$Cell_Index, "_S2", sep = "")

# binding both coldata from the 2 batches
coldata <- rbind(coldata1, coldata2)

SevenBridges Output

SevenBridges pre-processing pipeline report:

sb <- data.frame(Batch = c("S1", "S2"),
                 Number.of.Cells = c(dim(data1)[1], dim(data2)[1]),
                 Number.of.Features = c(dim(data1[, -1])[2], dim(data2[, -1])[2]),
                 Median.Number.of.Features = c(median(rowSums(data1 > 0)), median(rowSums(data2 > 0))),
                 stringsAsFactors = F)

DT::datatable(sb, 
              class = 'compact stripe hower',
              options = list(
                dom = "Bfrtip",
                scrollX = TRUE,
                paging = FALSE,
                searching = FALSE,
                info = FALSE,
                ordering = FALSE,
                columnDefs = list(list(className = 'dt-center', targets = 0:3))), rownames = FALSE)

Sample IDs

Number of cells in each sample:

table(coldata$Sample_Name)

## 
##   CARKOafter  CARKObefore     CARafter    CARbefore      KOafter     KObefore 
##         4554         2905         4831         3901         4836         3565 
##    Multiplet      NTafter     NTbefore Undetermined 
##         4128         4586         4052          221

We define a new column containing the group of samples:

• BeforeCoc = Before
• AfterCoc = After
• NA = Multiplets + Undetermined

coldata <- coldata %>%
  dplyr::mutate(group = ifelse(grepl("before", Sample_Name), "BeforeCoc",
    ifelse(grepl("after", Sample_Name), "AfterCoc", NA)))

p1 <- coldata %>%
  ggplot(aes(x = batch, fill = group)) + 
  geom_bar() +
  labs(x = "Batch", fill = "Group")

p2 <- coldata %>%
  ggplot(aes(x = batch, fill = Sample_Name)) + 
  geom_bar() +
  labs(x = "Batch", fill = "Sample Name")

p3 <- coldata %>%
  ggplot(aes(x = group, fill = Sample_Name)) + 
  geom_bar() +
  labs(x = "Group", fill = "Sample Name")

p1 + p2 + p3

coldata %>%
  ggplot(aes(x = batch, fill = group, color = Sample_Name)) +
  geom_bar() +
  labs(x = "Batch", fill = "Group", color = "Sample Name")

Pre-processing Workflow

Filtering

data_t <- data[, -1] %>% t()
rownames(data_t) <- colnames(data[, -1])
colnames(data_t) <- data$Cell_Index

sce <- SingleCellExperiment(assays = list(counts = data_t), colData = coldata)
rowData(sce)$Type <- ifelse(grepl("abo", rownames(sce), ignore.case = TRUE), "Antibody Capture", "Gene Expression") # 'Type' is either gene or antibody

# removing the cells with multiplet or undetermined assignment
sce <- sce[, !coldata$invalid]

coldata <- coldata %>%
  dplyr::filter(!invalid)

# splitting protein from RNA
sce_split <- splitAltExps(sce, rowData(sce)$Type)
# altExpNames(sce_split)

Overall Quality Metrics

Before Filtering

# ADT vs RNA counts
rna_numi <- colSums(counts(sce_split))
adt_numi <- colSums(counts(altExp(sce_split)))
nUMIdf <- data.frame(rna_nUMI = rna_numi,
                     adt_nUMI = adt_numi, 
                     Batch = sce_split$batch)

ggplot(nUMIdf, aes(x = log10(adt_nUMI+1), y = log10(rna_nUMI+1))) + 
  geom_point(alpha = 0.1, size = 1) +
  geom_density_2d(color = "orange") +
  facet_grid(.~Batch) +
  labs(x = "ADT Counts", y = "RNA Counts")

Computing QC metrics

outliers <- perCellQCMetrics(sce_split)
libsize.drop <- isOutlier(outliers$sum, log = TRUE, type = "both", nmads = 3)
feature.drop <- isOutlier(outliers$detected, log = TRUE, type = "lower", nmads = 3)
print(data.frame(ByLibSize = sum(libsize.drop),
                 ByFeature = sum(feature.drop),
                 Discard = sum(libsize.drop | feature.drop)))

##   ByLibSize ByFeature Discard
## 1       130       196     322

outliers.discard <- libsize.drop | feature.drop

#outliers.discard <- quickPerCellQC(outliers)
#colSums(as.matrix(outliers.discard))

Ab Counts

ab.discard <- isOutlier(outliers$`altexps_Antibody Capture_detected`,
                        log = TRUE,
                        type = "lower",
                        min_diff = 1)
summary(ab.discard)

##    Mode   FALSE 
## logical   33230

hist(outliers$`altexps_Antibody Capture_detected`, 
     col = 'grey',
     main = "", xlab = "Number of detected ADTs")

QC Plots

colData(sce_split) <- cbind(colData(sce_split), outliers) # adding the outliers info (sum, detected) to the sce object coldata
sce_split$discard <- outliers.discard | ab.discard # adding the cells to be discarded to the sce object coldata

plotColData(sce_split, x = "detected", y = "sum", colour_by = "discard", other_fields=c("batch")) +
  facet_wrap(~batch) +
  theme(panel.border = element_rect(color = "grey")) + 
  labs(x = "Number of Features Detected", y = "Library Size")

gridExtra::grid.arrange(
  plotColData(sce_split, x="batch", y="sum", colour_by="discard") +
    scale_y_log10() + ggtitle("Library Size"),
  plotColData(sce_split, x="batch", y="detected", colour_by="discard") +
    scale_y_log10() + ggtitle("Genes Detected") + labs(y = "genes detected"),
  plotColData(sce_split, x="batch", y="altexps_Antibody Capture_detected", colour_by="discard") +
    scale_y_log10() + ggtitle("Antibodies Detected") + labs(y = "antibodies detected"),
  ncol=1
)

gridExtra::grid.arrange(
  plotColData(sce_split, y="sum", colour_by="discard") +
    scale_y_log10() + ggtitle("Library Size"),
  plotColData(sce_split, y="detected", colour_by="discard") +
    scale_y_log10() + ggtitle("Genes Detected") + labs(y = "genes detected"),
  plotColData(sce_split, y="altexps_Antibody Capture_detected", colour_by="discard") +
    scale_y_log10() + ggtitle("Antibodies Detected") + labs(y = "antibodies detected"),
  ncol=1
)

# ADT vs RNA counts
rna_numi <- colSums(counts(sce_split))
adt_numi <- colSums(counts(altExp(sce_split)))
nUMIdf <- data.frame(rna_nUMI = rna_numi,
                     adt_nUMI = adt_numi, 
                     Batch = sce_split$batch,
                     Discard = sce_split$discard)

ggplot(nUMIdf, aes(x = log10(adt_nUMI+1), y = log10(rna_nUMI+1), color = Discard)) + 
  geom_point(alpha = 0.1, size = 1) +
  geom_density_2d(color = "black") +
  facet_grid(.~Batch) +
  labs(x = "ADT Counts", y = "RNA Counts") +
  scale_color_manual(values = c("#0072B2", "#C4961A"))

# filtering; removing outliers
sce_split <- sce_split[, !sce_split$discard]
sce_split

## class: SingleCellExperiment 
## dim: 431 32908 
## metadata(0):
## assays(1): counts
## rownames(431): ADA ADGRE1 ... ZBTB16 ZNF683
## rowData names(2): Type scDblFinder.selected
## colnames(32908): 620436_S1 864792_S1 ... 434005_S2 662384_S2
## colData names(17): Cell_Index Sample_Tag ... total discard
## reducedDimNames(0):
## mainExpName: Gene Expression
## altExpNames(1): Antibody Capture

table(coldata$Sample_Name)

## 
##  CARKOafter CARKObefore    CARafter   CARbefore     KOafter    KObefore 
##        4554        2905        4831        3901        4836        3565 
##     NTafter    NTbefore 
##        4586        4052

Normalization

# RNA

## size factors
lib.sf <- librarySizeFactors(sce_split)
summary(lib.sf)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.2508  0.6545  0.8914  1.0000  1.2232  3.6268

lib.sf.df <- data.frame(lib.sf, Group = colData(sce_split)$group)

## plotting size factors
hist(log10(lib.sf), xlab="Log10[Size Factor]", col='grey80')

ggplot(lib.sf.df, aes(log10(lib.sf))) +  
  geom_histogram() + 
  facet_wrap(~Group, scales = "free") + 
  labs(x = "Log10[Size Factor]")

## normalization
sce_split <- logNormCounts(sce_split)
#assayNames(sce_split)

# ADT

clr_norm <- function(x) {
  return(log1p(x = x / (exp(x = sum(log1p(x = x[x > 0]), na.rm = TRUE) / length(x = x)))))
}

logcounts(altExp(sce_split)) <- apply(counts(altExp(sce_split)), 2, clr_norm)

Variance Modeling

# determining variance components
dec <- modelGeneVar(sce_split)
dec[order(dec$bio, decreasing = TRUE), ]

## DataFrame with 431 rows and 6 columns
##             mean     total      tech       bio     p.value         FDR
##        <numeric> <numeric> <numeric> <numeric>   <numeric>   <numeric>
## CCL4     4.63735   3.88330  0.693030   3.19027 5.54825e-39 1.05417e-36
## IL32     5.01111   3.72651  0.661392   3.06512 1.75874e-39 6.68321e-37
## CCL3     3.39874   3.09598  0.954317   2.14166 1.13953e-10 1.44341e-08
## FCGR3A   2.58750   3.01710  1.070061   1.94704 1.36328e-07 1.03609e-05
## CCL5     3.69963   2.75220  0.883269   1.86893 1.12361e-09 1.06743e-07
## ...          ...       ...       ...       ...         ...         ...
## PTPRC    2.81000  0.637744   1.03353 -0.395788    0.860389    0.892215
## LYN      1.44121  0.739996   1.14443 -0.404430    0.840995    0.883346
## DOCK8    2.03512  0.741591   1.17331 -0.431720    0.850759    0.888155
## RPN2     2.25402  0.706351   1.14847 -0.442122    0.861652    0.892215
## LAMP1    2.41732  0.658019   1.11478 -0.456757    0.876517    0.902646

fit <- metadata(dec)
plot(fit$mean, fit$var, 
     xlab = "Mean of log-expression", ylab = "Variance of log-expression")
points(fit$mean, fit$var, col = "red", pch = 16)
curve(fit$trend(x), col = "dodgerblue", add = TRUE, lwd = 2)

Distribution of ADT Markers

Looking at the distribution of ADT markers (after normalization):

p <- do.call(plot_grid, c(lapply(rownames(logcounts(altExp(sce_split))), function(adt_name){
  ggplot(as.data.frame(logcounts(altExp(sce_split))[adt_name, ]), aes_string(x=logcounts(altExp(sce_split))[adt_name, ])) +
    geom_histogram(color="black",
                   fill = "black",
                   breaks=seq(0, 4.5, by=0.10)) + 
    xlab(adt_name) + 
    theme(axis.title.x = element_text(size = 8, face="bold"))}), ncol = 4))

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## i Please use tidy evaluation ideoms with `aes()`

print(p)

We define a new column containing the condition of samples:

• NT = NTbefore + NTafter
• KO = KObefore + KOafter
• CAR = CARbefore + CARafter
• CARKO = CARKObefore + CARKOafter

# creating new colData
sce_split$CAR <- ifelse(sce_split$Sample_Name %in% c("CARKOafter", "CARKObefore", "CARafter", "CARbefore"), "Yes", "No")
sce_split$KO <- ifelse(sce_split$Sample_Name %in% c("CARKOafter", "CARKObefore", "KOafter", "KObefore"), "Yes", "No")
sce_split$condition <- ifelse(sce_split$Sample_Name %in% c("CARKOafter", "CARKObefore"), "CARKO",
                                 ifelse(sce_split$Sample_Name %in% c("CARafter", "CARbefore"), "CAR",
                                        ifelse(sce_split$Sample_Name %in% c("KOafter", "KObefore"), "KO", "NT")))
sce_split$condition <- factor(sce_split$condition, levels = c("NT", "KO", "CAR", "CARKO"))

saveRDS(sce_split, "intermediate_files/D1/sce_donor1.Rds")

Processing - Downstream Analysis

Before

sce_split <- readRDS("intermediate_files/D1/sce_donor1.Rds")
sce_before <- sce_split[, colData(sce_split)$group == "BeforeCoc"]
sce_counts <- rbind(counts(sce_before), counts(altExp(sce_before)))
sce_logcounts <- rbind(logcounts(sce_before), logcounts(altExp(sce_before)))
sce <- SingleCellExperiment(assays = list(counts = sce_counts, logcounts = sce_logcounts))
rowData(sce)$Type <- ifelse(grepl("abo", rownames(sce), ignore.case = TRUE), "Antibody Capture", "Gene Expression")
sce$condition <- sce_before$condition
sce$condition <- gsub("NT", "NT-NK", sce$condition)
sce$condition <- gsub("CAR", "CAR33-NK", sce$condition)
sce$condition <- gsub("KO", "KLRC1.ko-NK", sce$condition)
sce$condition <- gsub("CAR33-NKKLRC1.ko-NK", "CAR33-KLRC1.ko-NK", sce$condition)
sce$condition <- factor(sce$condition, levels = c("NT-NK", "KLRC1.ko-NK", "CAR33-NK", "CAR33-KLRC1.ko-NK"))

set.seed(42L)

sce <- runPCA(sce, ncomponents=20)
sce <- runTSNE(sce, dimred="PCA")

snn_graph <- buildSNNGraph(sce, k=50, use.dimred='PCA')
igraph_clusters <- igraph::cluster_louvain(snn_graph)$membership
sce$igraph_lbls <- as.factor(igraph_clusters)

p1 <- plotReducedDim(sce, dimred = "TSNE", by_exprs_values = "logcounts", colour_by = "condition", other_fields = "condition") +
  facet_wrap(~condition, ncol = 4) + 
  scale_color_manual(values = c("#a0a0a3", "#90bff9", "#f2b77d", "#ff8080")) +
  theme_bw() + 
  theme(strip.text.x = element_text(size = 12), legend.title = element_blank(), legend.position="bottom", legend.text = element_text(size=10))

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

p2 <- plotReducedDim(sce, dimred = "TSNE", by_exprs_values = "logcounts", colour_by = "igraph_lbls", other_fields = "igraph_lbls") +
  facet_wrap(~igraph_lbls, ncol = 6) + 
  scale_color_manual(values = c("green", "orange", "purple", "navy", "cyan", "maroon")) +
  theme_bw() + 
  theme(strip.text.x = element_text(size = 12), legend.position="bottom", legend.text = element_text(size=10)) +
  guides(colour = guide_legend(nrow = 1)) +
  labs(color = "cluster")

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

plot_grid(p1 , p2, ncol = 1, align = "hv", labels = c("A", "B"))

## Warning: Graphs cannot be vertically aligned unless the axis parameter is set.
## Placing graphs unaligned.

ggsave("fig3a.tiff", p1, height = 2.5, width = 7.5)
ggsave("fig3c.tiff", p2, width = 7, height = 2)

ggsave("fig3a.pdf", p1, height = 2.5, width = 7.5)
ggsave("fig3c.pdf", p2, width = 7, height = 2)

p3 <- data.frame(colData(sce)) %>%
  group_by(igraph_lbls) %>%
  summarize(n = n()) %>%
  mutate(freq = n / sum(n) * 100) %>% 
  ggplot(aes(x = igraph_lbls, y = freq, fill = igraph_lbls)) +
  geom_col(position = "dodge", color = "black") + 
  scale_fill_manual(values = c("green", "orange", "purple", "navy", "cyan", "maroon"))+ 
  labs(x = "cluster", y = "cluster relative abundance (%)", fill = "cluster") +
  theme_bw() +
  theme(legend.position="bottom", legend.text = element_text(size=10)) +
  guides(fill = guide_legend(nrow = 1))

p4 <- data.frame(colData(sce)) %>%
  group_by(condition, igraph_lbls) %>%
  summarize(n = n()) %>%
  ggplot(aes(x = igraph_lbls, y = n, fill = condition)) +
  geom_col(position = "fill", color = "black") + 
  scale_fill_manual(values = c("gray", "blue", "orange", "red")) + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
  labs(title = "Cluster Abundance", x = "Cluster Number", y = "Cluster Relative Abundance")

## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.

p5 <- data.frame(colData(sce)) %>%
  group_by(condition, igraph_lbls) %>%
  summarize(n = n()) %>%
  mutate(freq = n / sum(n)) %>% 
  ggplot(aes(x = igraph_lbls, y = freq, fill = condition)) +
  geom_col(position = "dodge", color = "black") + 
  scale_fill_manual(values = c("#a0a0a3", "#90bff9", "#f2b77d", "#ff8080")) + 
  labs(x = "cluster", y = "condition relative abundance") +
  theme_bw() +
  theme(legend.position="bottom", legend.text = element_text(size=10)) +
  guides(fill = guide_legend(nrow = 1))

## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.

plot_grid(p3, p4, p5, ncol = 1, align = "hv", labels = c("A", "B", "C"))

## Warning: Graphs cannot be vertically aligned unless the axis parameter is set.
## Placing graphs unaligned.

## Warning: Graphs cannot be horizontally aligned unless the axis parameter is set.
## Placing graphs unaligned.

ggsave("fig3d.tiff", p3, height = 4, width = 4)
ggsave("fig3e.tiff", p5, height = 4, width=6)

ggsave("fig3d.pdf", p3, height = 4, width = 4)
ggsave("fig3e.pdf", p5, height = 4, width=6)

rownames(sce) <- gsub("TCR.alpha_beta.TRA_TRB.AHS0078.pAbO", "TCR_alpha_beta.TRA_TRB.AHS0078.pAbO", rownames(sce))
rownames(sce) <- gsub("TCR.gamma_delta.11F2.TRG_TRD.AHS0142.pAbO", "TCR_gamma_delta.11F2.TRG_TRD.AHS0142.pAbO", rownames(sce))
rownames(sce[grep("abo", rownames(sce), ignore.case = T, value = T), ]) <- gsub(".AH.*", "", rownames(sce[grep("abo", rownames(sce), ignore.case = T, value = T), ]))

Each dot represents a feature (gene or surface protein) within the respective pathway:

'%nin%' <- Negate('%in%')

markers <- c("HLA.DQA1", "HLA.DRA", "HLA.DQB1", "HLA.DPB1", "CD74",
             "TRAT1", "CD70", "CD86", "IL2RA", "CD25.2A3.IL2RA", "CD5", "CD5.UCHT2.CD5", "CD6",
             "CD45RA.HI100.PTPRC", "CD3D", "CD3E", "CD3G", "CD3.SK7.CD3E", "CD8A", "CD8B", "CD8.SK1.CD8A",
             "FCGR3A", "CD16.B73.1.FCGR3A_FCGR3B", "CD11b.ICRF44.ITGAM",
             "PDCD1", "CD279.EH12.1.PDCD1", "CTLA4", "TIGIT", "Tim3.HAVCR2", "CD94.KLRD1", "ZNF683", "CD38", "CD38.HIT2.CD38")

score <- data.frame(matrix(NA, nrow = length(levels(sce$igraph_lbls)), ncol = length(markers)))
rownames(score) <- levels(sce$igraph_lbls)
colnames(score) <- markers

for(m in markers){
  for(c in levels(sce$igraph_lbls)){
    mean.all <- counts(sce)[m, sce$igraph_lbls %nin% c] %>% mean()
    mean.cluster <- counts(sce)[m, sce$igraph_lbls %in% c] %>% mean()
    if(mean.cluster == 0){
      score[c, m] <- -1
    } else {
      score[c, m] <- log10(mean.cluster/mean.all)
    }
  }
}

score <- score %>% rownames_to_column('cluster') %>% pivot_longer(cols = -cluster)
colnames(score) <- c("cluster", "markers", "value")
score$pathway <- case_when(
  score$markers %in% c("HLA.DQA1", "HLA.DRA", "HLA.DQB1", "HLA.DPB1", "CD74") ~ "Downstream TCR Signaling",
  score$markers %in% c("TRAT1", "CD70", "CD86", "IL2RA", "CD25.2A3.IL2RA", "CD5", "CD5.UCHT2.CD5", "CD6") ~ "Lymphocyte Activation",
  score$markers %in% c("CD45RA.HI100.PTPRC", "CD3D", "CD3E", "CD3G", "CD3.SK7.CD3E", "CD8A", "CD8B", "CD8.SK1.CD8A") ~ "Naive Lymphocyte Markers",
  score$markers %in% c("FCGR3A", "CD16.B73.1.FCGR3A_FCGR3B", "CD11b.ICRF44.ITGAM") ~ "Mature NK Cells",
  score$markers %in% c("PDCD1", "CD279.EH12.1.PDCD1", "CTLA4", "TIGIT", "Tim3.HAVCR2", "CD94.KLRD1", "ZNF683", "CD38", "CD38.HIT2.CD38") ~ "Checkpoints & Suppressive Regulators"
  )

p10 <- ggplot(score, aes(x = cluster, y = value, fill = cluster)) +
  facet_wrap(~pathway, ncol = 5) +
  geom_dotplot(binaxis = 'y', stackdir = 'center', stackratio = 1, dotsize = 2, binwidth = 0.05) +
  geom_hline(yintercept = 0,  col = "red", size = 0.5) + 
  stat_summary(fun = median, geom = "point", shape = 18, size = 1.65, color = "black") +
  scale_fill_manual(values = c("green", "orange", "purple", "navy", "cyan", "maroon")) +
  theme_bw()+ 
  theme(legend.position="none", strip.text.x = element_text(size = 10)) + 
  ylab("differential expression score")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## i Please use `linewidth` instead.

p10

ggsave("fig3f.tiff", p10, width = 14, height = 3)
ggsave("fig3f.pdf", p10, width = 14, height = 3)

sce_split <- splitAltExps(sce, rowData(sce)$Type)
rownames(altExp(sce_split)) <- gsub("\\..*", "", rownames(altExp(sce_split)))
altExp(sce_split)$igraph_lbls <- sce_split$igraph_lbls

supp <- plotDots(altExp(sce_split), features = rownames(altExp(sce_split)), group = "igraph_lbls",
                 color = c("darkblue", "lightblue", "gray", "yellow", "orange", "darkorange", "red", "darkred")) +
  theme(legend.position="bottom", legend.box="vertical", legend.margin=margin())

ggsave("fig3supp.tiff", supp, width = 6, height = 12)
ggsave("fig3supp.pdf", supp, width = 6, height = 12)

Session

date()

## [1] "Thu Jan 26 10:46:21 2023"

sessionInfo()

## R version 4.2.2 Patched (2022-11-10 r83330)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.5 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=C                 LC_NUMERIC=C              
##  [3] LC_TIME=de_DE.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=C             
##  [7] LC_PAPER=de_DE.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] ggrepel_0.9.2               edgeR_3.38.4               
##  [3] limma_3.52.4                pheatmap_1.0.12            
##  [5] cowplot_1.1.1               scDblFinder_1.10.0         
##  [7] patchwork_1.1.2             scater_1.24.0              
##  [9] scran_1.24.1                scuttle_1.6.3              
## [11] SingleCellExperiment_1.18.1 SummarizedExperiment_1.26.1
## [13] Biobase_2.56.0              GenomicRanges_1.48.0       
## [15] GenomeInfoDb_1.32.4         IRanges_2.30.1             
## [17] S4Vectors_0.34.0            BiocGenerics_0.42.0        
## [19] MatrixGenerics_1.8.1        matrixStats_0.63.0         
## [21] forcats_0.5.2               stringr_1.5.0              
## [23] dplyr_1.0.10                purrr_1.0.1                
## [25] readr_2.1.3                 tidyr_1.3.0                
## [27] tibble_3.1.8                ggplot2_3.4.0              
## [29] tidyverse_1.3.2            
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.4.1              backports_1.4.1          
##   [3] systemfonts_1.0.4         igraph_1.3.5             
##   [5] BiocParallel_1.30.4       crosstalk_1.2.0          
##   [7] digest_0.6.31             htmltools_0.5.4          
##   [9] viridis_0.6.2             fansi_1.0.4              
##  [11] magrittr_2.0.3            ScaledMatrix_1.4.1       
##  [13] googlesheets4_1.0.1       cluster_2.1.4            
##  [15] tzdb_0.3.0                Biostrings_2.64.1        
##  [17] modelr_0.1.10             timechange_0.2.0         
##  [19] colorspace_2.1-0          rvest_1.0.3              
##  [21] textshaping_0.3.6         haven_2.5.1              
##  [23] xfun_0.36                 crayon_1.5.2             
##  [25] RCurl_1.98-1.9            jsonlite_1.8.4           
##  [27] glue_1.6.2                gtable_0.3.1             
##  [29] gargle_1.2.1              zlibbioc_1.42.0          
##  [31] XVector_0.36.0            DelayedArray_0.22.0      
##  [33] BiocSingular_1.12.0       scales_1.2.1             
##  [35] DBI_1.1.3                 Rcpp_1.0.10              
##  [37] isoband_0.2.7             viridisLite_0.4.1        
##  [39] dqrng_0.3.0               rsvd_1.0.5               
##  [41] DT_0.27                   metapod_1.4.0            
##  [43] htmlwidgets_1.6.1         httr_1.4.4               
##  [45] RColorBrewer_1.1-3        ellipsis_0.3.2           
##  [47] pkgconfig_2.0.3           XML_3.99-0.13            
##  [49] farver_2.1.1              sass_0.4.5               
##  [51] dbplyr_2.3.0              locfit_1.5-9.7           
##  [53] utf8_1.2.2                tidyselect_1.2.0         
##  [55] labeling_0.4.2            rlang_1.0.6              
##  [57] munsell_0.5.0             cellranger_1.1.0         
##  [59] tools_4.2.2               cachem_1.0.6             
##  [61] xgboost_1.7.3.1           cli_3.6.0                
##  [63] generics_0.1.3            broom_1.0.2              
##  [65] evaluate_0.20             fastmap_1.1.0            
##  [67] ragg_1.2.5                yaml_2.3.7               
##  [69] knitr_1.41                fs_1.6.0                 
##  [71] sparseMatrixStats_1.8.0   xml2_1.3.3               
##  [73] compiler_4.2.2            rstudioapi_0.14          
##  [75] beeswarm_0.4.0            reprex_2.0.2             
##  [77] statmod_1.5.0             bslib_0.4.2              
##  [79] stringi_1.7.12            highr_0.10               
##  [81] lattice_0.20-45           bluster_1.6.0            
##  [83] Matrix_1.5-3              vctrs_0.5.2              
##  [85] pillar_1.8.1              lifecycle_1.0.3          
##  [87] jquerylib_0.1.4           BiocNeighbors_1.14.0     
##  [89] data.table_1.14.6         bitops_1.0-7             
##  [91] irlba_2.3.5.1             rtracklayer_1.56.1       
##  [93] R6_2.5.1                  BiocIO_1.6.0             
##  [95] gridExtra_2.3             vipor_0.4.5              
##  [97] codetools_0.2-18          MASS_7.3-58.1            
##  [99] assertthat_0.2.1          rjson_0.2.21             
## [101] withr_2.5.0               GenomicAlignments_1.32.1 
## [103] Rsamtools_2.12.0          GenomeInfoDbData_1.2.8   
## [105] parallel_4.2.2            hms_1.1.2                
## [107] grid_4.2.2                beachmat_2.12.0          
## [109] rmarkdown_2.20            DelayedMatrixStats_1.18.2
## [111] googledrive_2.0.0         Rtsne_0.16               
## [113] lubridate_1.9.1           ggbeeswarm_0.7.1         
## [115] restfulr_0.0.15

knitr::knit_exit()